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Jia  Zhao*and  Qi  Wang* 


Abstract 

Biofilms  are  known  to  be  more  persistent  to  antimicrobial  treatment  than  planktonic  bacterial 
cells.  In  addition  to  the  protective  extracellular  matrix  formed  primarily  by  expopoly saccharides 
(EPS),  one  of  the  current  point  of  views  on  this  issue  is  that  there  may  exist  a  small  portion  of  phe¬ 
notypic  bacterial  variants,  known  as  persisters,  which  are  invulnerable  to  antimicrobial  agents,  while 
the  majority  of  the  bacterial  cells  is  susceptible  to  antimicrobial  agents.  In  this  paper,  a  3D  hydrody¬ 
namic  model  for  spatially  heterogeneous  biofilms  based  on  the  phase  field  formulation  is  proposed 
and  applied  to  analyze  the  mechanism  of  antimicrobial  persistence  of  biofilms  by  acknowledging  the 
existence  of  persisters  and  susceptible  cells  in  the  total  population  of  bacteria.  A  numerical  scheme 
is  devised  to  solve  the  model  consisting  of  partial  differential  equations,  which  is  implemented  on 
graphic  processing  units  (GPUs)  for  high  performance  computing,  in  3-D  space  and  time.  Antimicro¬ 
bial  treatment  in  an  infinitely  long  quiescent  water  channel  and  in  a  short  water  tube  under  inflow  and 
outflow  boundary  conditions  are  simulated  using  the  new  model,  in  which  multiple  dosing  locations 
and  strategies  are  investigated.  The  model  demonstrates  the  internal  spatial  and  temporal  distribution 
of  bacteria,  EPS,  nutrient  and  antimicrobial  agents,  providing  a  useful  tool  for  analyzing  the  mecha¬ 
nism  of  biofilm  persistence  to  antimicrobial  agents  in  an  aqueous  environment.  The  numerical  result 
also  confirms  that  the  periodic  dosing  strategy  is  more  effective  than  the  constant  dosing  strategy  in 
disinfecting  biofioms. 


1  Introduction 

In  nature,  as  long  as  bacteria  colonize  on  moisture  surfaces,  biofilms  will  likely  form,  which  are  con¬ 
sisted  of  the  micro-organisms  aggregated  by  bacteria  colonies  along  with  their  self-produced,  glue-like 
polysacharride  matrix,  known  as  the  extracellular  polymeric  substance  (EPS).  It’s  commonly  perceived 
by  the  medical  community  that  biofilms  are  responsible  for  many  diseases  or  ailments  associated  with 
chronic  infections,  evidenced  by  the  survey  that  biofilms  are  present  on  the  removed  tissue  of  80%  of  pa¬ 
tients  undergoing  surgery  for  chronic  sinusitis  [36].  Unlike  a  planktonic  bacterium,  biofilms  are  always 
hard  to  be  eradicated  by  the  standard  antimicrobial  treatment  [28],  which  perhaps  explains  the  frequent 
relapse  of  chronic  diseases  or  ailments  associated  with  biofilms. 

Thus,  an  understanding  of  the  mechanism  that  underlies  biofilm  persistence  to  antimicrobial  agents 
can  greatly  enhance  therapeutic  treatment  of  diseases  related  with  biofilms.  Intensive  research  efforts 
have  been  carried  out,  primarily  experimental,  to  try  to  understand  biofilm  dynamics,  yet  little  is  fully 
known.  Readers  may  refer  to  the  review  papers  [12]  and  [28]  for  overviews  of  current  advances  in  the 
treatment  of  biofilms.  What  is  known  to  us  via  experimental  evidence  is  that  not  only  one,  but  many 
factors  can  contribute  to  the  biofilm-antimicrobial  agent  interaction.  Among  them,  one  essential  factor 
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is  the  existence  of  persister  cells  within  the  biofilm  colony,  which  are  consisted  of  a  small  portion  of 
dormant  bacterial  variants  and  are  highly  tolerant  to  antimicrobial  agents  [3].  Contrasting  to  persister 
cells,  the  other  bacteria  are  collectively  called  susceptible  bacteria. 

From  the  clinical  point  of  view,  understanding  the  mechanism  of  persister  formation  would  be  es¬ 
sential  for  biofilm  control  and  thereby  impact  on  the  treatment  of  diseases  responsible  by  biofilms.  For 
review  papers  on  mechanisms  underlying  the  persister  formation,  readers  are  referred  to  the  works  by 
Kim  Lewis  [28]  and  [27].  As  dormant  variants  of  regular  bacterial  cells,  which  doesn’t  undergo  genetic 
changes,  it  is  convinced  that  persisters  are  converted  from  regular  cells  due  to  stresses  [3],  such  as  nutri¬ 
ent  depletion  [5],  hydrodynamical  shear,  existence  of  antimicrobial  agents  [33]  and  so  on.  Later,  when 
the  environment  is  tolerable,  say  nutrient  is  sufficient  or  the  concentration  of  antimicrobial  agents  drops 
under  a  certain  threshold  value,  biofilms  can  relapse  [6],  which  implies  that  persisters  convert  back  into 
susceptible  bacteria  for  regrowth. 

Factoring  in  persister  formation,  people  have  conducted  research  on  therapeutic  treatment  of  dis¬ 
eases  induced  by  biofilms.  The  review  paper  [39]  provides  some  information  about  control  strategies  for 
biofilms.  Concerning  dosing  strategies  of  antimicrobial  agents,  there  exists  an  evidence  that  a  concen¬ 
trated  dose  of  biocide  is  more  effective  than  using  a  prolonged  dose  of  a  lower  concentration  [20].  In 
addition,  dosing  by  shocks  is  more  effective  than  dosing  in  a  persistent  manner  [19].  But  to  the  best  of  our 
knowledge,  there  does  not  exist  any  optimal  strategies  derived  for  biofilm  control  or  disease  treatment 
caused  by  biofilms  so  far.  Currently,  the  environmental  impact  of  biocide  or  side-effect  of  antibiotics 
have  become  common  concerns,  which  makes  the  derivation  of  an  optimal  strategy  even  harder. 

Besides  the  formation  of  persister  bacteria,  extra  cellular  substances  (EPS)  also  play  a  role  in  the 
failure  of  antimicrobial  treatment  of  biofilms.  Basically,  EPS  is  believed  to  act  as  a  protective  barrier 
to  prevent  the  antimicrobial  agents  from  penetrating  deep  into  the  biofilm  region,  either  by  reacting 
with  antimicrobial  agents  [38]  or  by  simply  slowing  down  the  diffusion  rate  via  its  densely  distributed 
network  meshes  [41].  Readers  can  find  more  details  about  the  role  of  EPS  in  biofilm  structure  and 
function  from  [17]  and  [16]. 

From  the  mathematical  perspective,  several  models  have  been  proposed  trying  to  interpret  experi¬ 
mental  observations  on  biofilm  structures  and  function.  For  a  review  of  mathematical  models  of  biofilms, 
interested  readers  are  referred  to  [25]  and  [44]  for  more  details.  Mathematical  models  at  different  length 
scales  have  also  been  devised,  i.e.,  microscopic  scale  models  (agent  based  models)  [1],  mesoscopic  scale 
models  (kinetic  theory  models)  [45]  and  macroscopic  scale  models  (continuum  theories)  [47],  based  on 
the  specific  issue  that  one  was  interested  in.  Recently,  modeling  biofilms  as  multiphase  complex  fluids 
has  emerged  as  a  promising  approach  to  address  some  complex  and  intriguing  issues  associated  with 
biofilm  dynamics  [46,47],  where  bacteria  are  regarded  as  colloids  and  the  biofilm  as  a  polymer  gel.  In 
such  an  approach,  a  complex  fluid  model  can  be  devised  to  analyze  the  structure  formation  and  function 
of  biofilms  in  a  hydrodynamical  setting.  In  this  research  direction,  the  work  of  Wang  et  al.  represents 
some  latest  development  [47,48]. 

On  the  issue  of  biofilm  persistence  to  antimicrobial  agents,  simple  mathematical  models  have  also 
been  developed  to  test  certain  mechanisms  for  persister  formation  based  on  the  experimental  evidence 
that  supports  the  concept  of  persisters.  For  instance  in  [35],  the  author  claims  using  a  simple  mathe¬ 
matical  model  that  persister  formation  can  lead  to  higher  bacterial  persistence  to  antimicrobial  agents 
than  those  grown  in  plancktonic  culture.  In  [23],  a  3D  agent-based  model  for  biofilm  dynamics  under 
antimicrobial  treatment  was  developed,  in  which  it  showed  that  substrate  limitation  can  contribute  to  the 
persistence  of  biofilms  to  antimicrobial  agents.  Cogan  has  worked  on  possible  mechanisms  of  persister 
formation  using  time-dependent,  but  spatially  homogeneous  models  recently  [8],  [11],  [26]. 

Models  on  dosing  strategies  for  treating  diseases  caused  by  biofilms  have  also  been  proposed.  Cogan 
discussed  effective  dosing  strategies  using  simple  mathematical  models  in  [8,26].  In  [10],  he  discussed 
the  effect  of  periodic  disinfection  by  a  one-dimensional  mathematical  model.  In  [43],  the  adaptive  re- 
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sponse  to  dosing  protocols  for  biofilm  control  has  also  been  analyzed,  which  provides  some  sufficient 
conditions  for  eradicating  biofilms  under  the  constant  dosing  approach.  In  addition,  models  analyzing 
other  impact  factors,  which  may  contribute  to  biofilm’s  persistence  to  antimicrobial  agents,  have  also 
been  proposed.  For  instance,  the  author  in  [13]  tried  to  analyze  and  simulate  diffusive  resistance  of 
bacterial  biofilms  to  penetration  of  antibiotics. 

To  our  best  knowledge,  there  has  yet  been  a  mathematical  model  which  takes  into  account  both  the 
hydrodynamic  effect  and  the  spatio-temporal  heterogeneous  structures  of  biofilms  in  full  three  dimen¬ 
sion  in  space  and  time.  However,  these  factors  are  important  in  biofilm  structure  formation  and  function, 
especially,  concerning  with  the  biofilm’s  persistence  to  antimicrobial  agents.  In  this  paper,  We  develop 
a  full  3-D  hydrodynamic  model  for  biofilms  of  multiple  bacterial  phenotypes;  in  particular,  we  limit  the 
phenotypes  to  the  persister  and  susceptible  type.  This  model  extends  our  previous  model  for  biofilm- 
solvent  mixture  [47]  by  distinguishing  between  the  persister  cell  and  the  susceptible  cell  when  biofilms 
are  treated  with  antimicrobial  agents.  In  this  model,  the  interplay  among  the  various  biomass  components 
such  as  various  bacterial  types,  EPS  and  solvent  is  carefully  accounted  for  both  hydrodynamically  and 
chemically  so  that  dead  bacteria  can  be  deteriorate  into  solvent  and  decomposed  into  EPS  simultaneously 
some  time  after  their  death.  Our  model  shows  that  the  dynamical  interaction  between  these  two  general 
phenotypes  can  impact  dramatically  on  the  overall  dynamics  of  the  biofilm.  The  most  distinctive  at¬ 
tribute  of  this  model  is  that  it  provides  the  spatio-temporal  resolution  that  can  resolve  more  details  about 
antimicrobial  action  against  biofilm  colonies  in  space  and  time  than  the  previous  models  can,  providing 
a  much  needed  insight  into  dynamics  of  antimicrobial  treatment. 

The  rest  of  this  paper  is  organized  as  follows:  In  section  two,  we  provide  the  formulation  of  the 
hydrodynamic  theory  for  the  biofilm  system  based  on  the  phase  field  formulation.  Then,  an  efficient  nu¬ 
merical  solver  for  this  coupled  fluid  matrix-biofilm  partial  differential  equation  system  by  semi-implicit 
finite  difference  strategy  is  proposed  in  section  three.  In  section  four,  numerical  simulations  and  results 
are  presented  and  discussed.  Finally,  we  summarize  the  result  and  draw  a  conclusion. 

2  Mathematical  Model  Formulation 

2.1  Notations 

We  model  the  biofilm  together  with  its  surrounding  aqueous  environment  as  a  mixture  of  complex 
fluids.  The  biofilm  is  consisted  of  the  biomass  immersed  in  solvent;  whereas  the  biomass  is  made  up 
of  bacteria  and  their  products  like  exopolysaccharide  (EPS).  Nutrient  and  antimicrobial  agents  are  small 
molecule  substances  dissolved  in  solvent.  Let  q !>&s  be  the  volume  fraction  of  the  bacteria  that  are  suscep¬ 
tible  to  antimicrobial  agents,  fi\yp  the  volume  fraction  of  the  bacteria  that  are  persistent  to  antimicrobial 
agents,  <fibd  the  volume  fraction  of  the  dead  bacteria,  and  cj)p  the  volume  fraction  of  EPS.  We  note  that 
the  EPS  is  the  product  of  the  live  bacteria  at  the  presence  of  nutrient.  Dead  bacteria  do  not  produce  any 
EPS  and,  over  times,  dead  bacteria  can  dissolve  into  solvent  and  part  of  them  will  be  released  as  EPS 
in  the  biofilm  colony  because  some  EPSs  are  attached  to  the  cell  membrane  at  the  moment  when  the 
cell  dies.  We  in  addition  denote  the  concentration  of  the  nutrient  and  the  antimicrobial  agent  as  c  and 
d ,  respectively.  To  make  contact  with  our  previous  binary  mixture  model  for  biofilms,  we  define  <fin  the 
volume  fraction  of  the  biomass,  which  consists  of  all  the  volume  fractions  for  the  bacteria  fib  as  well  as 
EPS  fip,  i.e., 

fin  =  fib  +  fip,  fib  ~  fibs  +  fibp  +  fibd- 

In  addition  to  the  volume  fractions  introduced  above,  the  volume  fraction  of  the  solvent  is  denoted  as 
fis.  Due  to  the  small  molecular  weigh  in  the  nutrient  and  antimicrobial  molecules,  we  will  not  explicitly 
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count  their  mass  and  volume  in  this  model.  The  incompressibility  of  the  complex  fluid  mixture  then 
implies  (j)s  +  c/)n  =  1. 

In  this  model,  we  make  a  simplifying  assumption  that  all  component  in  the  biomass  including  the 
bacteria  and  EPS  share  the  same  mass  density.  We  denote  pn  and  ps  the  density  of  biomass  and  solvent, 
and  vn  and  vs  the  velocity  of  biomass  and  solvent,  respectively.  Then,  the  volume  averaged  velocity  and 
density  are  given  respectively  by 

v  =  0nvn  T"  P  —  (frnPn  T  4>sPs- 

We  assume  the  bacteria,  regardless  whether  they  are  live  or  dead,  persisters  or  susceptibles,  and  the 
EPS  mix  with  the  solvent  owing  to  the  osmotic  pressure.  Then,  we  adopt  the  modified  mixing  free  energy 
introduced  in  (  [47],  [48])  and  denote  the  energy  density  by  f: 

/  =  y  fcsT||V0n||2  +  i2kBT  ln0„  +  (1  -  </>n)  ln(l  -  <f>n)  +  xM1  ~  4>n))  •  (1) 

This  is  the  modified  Flory-Huggins  mixing  free  energy  with  a  conformational  entropy  added,  in  which 
71,  72  parametrize  the  strength  of  the  conformational  entropy  and  bulk  mixing  free  energy,  respectively, 
X  is  the  mixing  parameter,  N  is  the  extended  polymerization  index  for  the  biomass,  ks  is  the  Boltzmann 
constant  and  T  is  the  absolute  temperature. 


2.2  Transport  equations  for  biomass  components 

Given  the  mixing  free  energy  density  f  in  (1),  the  ’’extended”  chemical  potentials  with  respect  to  each 
component  in  the  biomass  are  summarized  as  follows 

...  Sf  ..  Sf  Sf  _  Sf 

Pbs  r  /  5  Pbp  r  I  5  PJbd  r  /  5  Pp  r  /  * 

OCpbs  Otybp  Otftbd  S(pp 


The  transport  equation  for  the  volume  fraction  of  each  biomass  component  is  governed  by  a  reactive 
Cahn-Hilliard  equation, 


mfos  +  v  •  (<j)hs v) 

M&bp  +  V  •  (06pV) 
ffifod  +  V  •  (<?Wv) 
M0P  +  V-  (0pv) 


^  '  (^bsfibs^ Pbs)  T  pbs  5 

V  •  (XbpfibpV Pbp)  9bp-> 

V  •  (A bd<t>bdy Pbd)  +  9bd , 

V  •  (A p4>pS7 Pp)  +  9pi 


(2) 


where  9bs^9bp^9bd  and  gp  are  the  reactive  terms  (rates)  for  the  susceptible  bacteria,  persisters,  dead 
bacteria,  and  the  EPS  component,  respectively,  which  will  be  given  in  details  in  the  next  subsection, 
A bs,  A bp,  A bd  and  Xp  are  the  motility  parameters  for  the  transport  of  the  four  components,  which  can  be 
functions  of  the  biomass  volume  fractions. 


2.3  Transport  equations  for  nutrient  and  antimicrobial  agents 

In  this  model,  the  nutrient  is  treated  as  a  phantom  material,  in  which  its  mass  is  completely  neglected, 
but  the  chemical  effects  are  retained.  This  is  because  nutrient  is  consisted  of  small  molecule  materials 
compared  to  the  biomass,  which  are  made  up  of  the  bacterial  cells  and  large  molecular  EPS.  Specifically, 
oxygen  is  one  of  the  key  ingredients  in  nutrient  for  the  biofilm  formation. 
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The  governing  equation  for  the  nutrient  concentration  (c)  is  given  by  a  reaction-convection-diffusion 
equation  with  a  varying  diffusion  coefficient  and  reactive  term,  namely, 


d(4>sc) 

dt 


+  V  •  (cvs0s) 


V  •  (Ds(/>sVc)  +  gc, 


(3) 


where  Ds  is  the  diffusion  coefficient  and  gc  is  the  nutrient  consumption  term.  It  was  suggested  in  [40] 
that 

Dc, polymer  n 


where  DC:POiyrner  is  the  diffusion  rate  of  nutrient  in  a  fully  developed  biofilm  and  Dc  is  the  diffusion  rate 
in  the  pure  solvent.  In  particular,  Dc  is  measured  as  2  x  10-9  m-2s-1  at  temperature  25  °C  for  oxygen. 
Since  the  molecular  mass  of  oxygen  is  relatively  small,  we  propose  that  it  can  penetrate  the  EPS  and  the 
membrane  of  cells  equally.  Therefore,  the  diffusion  coefficient  of  oxygen  in  biofilms  can  be  formulated 
as 


Ds  =  Dc 


2(1  -  M 

2  +  (t>n 


where  the  second  term  on  the  right  represents  the  reduction  of  the  diffusion  rate  due  to  the  presence  of 
the  biomass. 

Like  the  nutrient,  an  antimicrobial  agent  can  also  be  modeled  as  a  phantom  material.  The  transport 
equation  of  it  is  proposed  as  follows 


+  V  •  (dvs(f>s)  =  V  •  ( De<t>s\7d )  +  gd,  (4) 

where  De  is  the  diffusion  coefficient  of  the  antimicrobial  agents  and  gd  is  the  reactive  term  including  the 
antimicrobial  injection  rate.  For  some  antimicrobial  agents,  such  as  the  Penicillin,  which  is  of  a  larger 
molecule,  the  biomass  (live  bacteria  and  EPS,  as  well  as  the  dead  cells)  can  prevent  the  antimicrobial 
agents  from  penetrating  deeper  into  the  biofilm  over  times.  To  quantify  this  diffusive  effect,  a  compara¬ 
tive  study  of  three  cases,  by  adjusting  the  diffusion  rate  of  antimicrobial  agents,  will  be  conducted  in  this 
paper.  The  simplest  case  is  proposed  as  follows 


Dei  =  Dd , 


(5) 


which  represents  an  isotropic  diffusion.  By  considering  the  fact  that  the  network  of  EPS  and  the  presence 
of  bacterial  cells  can  slow  down  the  penetration,  the  second  case  is  proposed  as  follows 


De  2  =  Dd 


2(1  ~  h) 
2  +  06 


05 


05  + 


(6) 


where  Dpr  is  a  parameter  in  the  Hinson  model  [22]  fitted  experimentally.  In  the  last  case,  we  assume 
that  both  EPS  and  various  bacteria  in  the  biofilm  can  slow  down  the  diffusion  of  antimicrobial  agents  at 
various  rates.  The  diffusion  rate  is  formulated  as  follows 


~j~~\  _  2 [1  (065  +  <fibp)\  05 

^  2  +  06^  +  06,  0S  +  ^* 

Upr 


(7) 


On  the  right  hand  side  of  the  formula  for  De 2  and  De 3,  the  second  term  represents  the  reduction  due  to  the 
physical  presence  of  bacteria,  and  the  third  term  is  an  empirical  fitting,  accounting  for  the  reduction  effect 
due  to  the  presence  of  EPS  alone  or  EPS  and  dead  bacteria  combined.  We  use  Dpr  —  0.02  throughout 
the  paper. 
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2.4  Hydrodynamic  equations  for  the  complex  fluid  mixture 


To  complete  this  biological  model,  the  balance  equation  for  the  averaged  velocity  and  for  the  mass 
density  is  given  below,  respectively.  For  the  average  velocity  v,  it’s  assumed  to  be  solenoidal.  Then,  the 
continuity  and  the  momentum  balance  equations  are  given,  respectively,  by 

p(<9tv  +  V  •  Vv)  =  V  •  ((f>bTb  +  (f)pTp  +  <j)sTs)  -  [Vp  +  7ifcBTV  •  (S7(f)n  0  V0„)], 


where  p  is  the  hydrostatic  pressure,  r&,  rp  and  rs  are  the  stress  tensor  due  to  bacteria,  the  EPS,  and 
solvent,  respectively.  The  last  term  in  the  momentum  equation  is  due  to  the  spatial  inhomogeneity  of  the 
biomass  distribution,  which  is  derived  from  the  least  action  principle  [47]. 

In  this  paper,  we  are  interested  in  the  growth  dynamics  of  the  biofilm,  whose  time  scale  is  significantly 
larger  than  the  relaxation  time  scale  in  the  EPS.  We  therefore  model  all  effective  biofilm  components  as 
viscous  fluids: 

Tb  —  2775D77,,  Tp  —  2fjpY)rn  ts  — 

where  775,  r]p  and  ps  are  the  viscosity  of  bacteria,  the  EPS  and  solvent,  respectively,  which  can  be  func¬ 
tions  of  volume  fractions  of  the  effective  biofilm  components.  Here,  the  rate  of  strain  tensors  for  each 
component  is  defined  respectively  by 

Dn  =  1  (Vv„  +  Vv^) ,  Ds  =  I  (Vvs  +  Vvf )  , 
where  we  assume  the  rates  of  strain  tensor  for  the  effective  biomass  components  are  the  same  given  by 


2.5  Reactive  kinetics  of  biomass  components  in  biofilms 

In  order  to  complete  this  model,  we  need  to  propose  reactive  kinetics  for  the  effective  biomass  com¬ 
ponents,  for  which  it  is  necessary  for  us  to  sort  out  the  relations  among  the  various  components  with 
respect  to  their  reactive  kinetics. 


2.5.1  Reactive  kinetics  of  susceptible  bacteria  and  persisters 


The  bacterial  growth  depends  on  the  nutrient  and  the  concentration  of  antimicrobial  agents  present  in 
the  biofilm  colony.  The  two  different  bacterial  phenotypes  (susceptible  and  persister  cells)  have  different 
growth  mechanisms;  thus  we  assume  that  both  the  susceptible  and  the  persister  cell  grow  on  their  own  at 
their  own  rates,  respectively.  For  both  bacterial  phenotypes,  death  rates  due  to  natural  causes  are  taken 
into  account  as  well.  In  addition,  both  susceptible  bacteria  and  persisters  can  be  killed  by  antimicrobial 
agents,  even  though  persisters  are  killed  in  a  much  slower  rate,  due  to  its  antimicrobial  persistence.  We 
allow  all  these  features  in  the  model  so  that  they  can  be  turned  on  and  off  depending  on  the  time  scale 
that  we  try  to  resolve.  It  is  perceived  that  susceptible  bacteria  and  persisters  can  be  converted  mutually, 
based  on  the  stage  of  their  growth  [28]  and  the  surrounding  environment  [8],  such  as  accessibility  to 
nutrient  and  antimicrobial  agents  etc. 

The  reactive  kinetics  for  these  two  types  of  live  bacteria  are  proposed,  respectively,  as  follows: 


9bs  ~  Ki+c( 1 


Vbs  \ _ ns _ U  rf)i  —1—  /)  rhu  (  rhsKsd  I  C’3 d  \1 


Vps^bp  V  K2sd+c 2  “r  K3+d  > 


ybp  —  K2+c 


(1  (t>bv^^^bP  +  ^sP^bs  bpS(j)bp  (Knd+c  +  d+K3 


ld+c  1  d+KsWOpi 


(9) 
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where  C2  and  C4  are  the  growth  rate  for  the  susceptible  and  the  persister,  respectively;  max?  4>bp, max 
are  the  maximum  volume  fraction  for  susceptible  bacteria  and  persisters,  respectively;  bsp  and  bps  are 
the  conversion  rates  between  the  two  types  of  bacteria;  r&s,  r&p  are  the  natural  death  rates  due  to  nutrient 
depletion  for  susceptible  bacteria  and  persisters,  respectively,  and  C3,  C12  are  the  death  rates  for  sus¬ 
ceptible  bacteria  and  persisters  due  to  antimicrobial  agents,  respectively.  Ki,K2,  Ksd ,  Kpci  are  half 

saturation  constants  in  the  Monod  models  adopted.  We  assume  r ^  and  C3  >>  C12  in  the  above 

model.  We  note  the  extra  term  ^  ^  in  g ^  in  the  susceptible  bacterial  growth  is  meant  to  rule  out  the 
effect  of  planktonic  susceptible  bacteria  in  biofilms. 

In  the  following,  we  propose  systematically  a  formulation  for  the  switch  function  bsp  and  bps  between 
the  two  types  of  bacteria,  based  on  the  experimental  description,  as  well  as  literature  on  mathematical 
interpretation.  It  is  basically  agreed  that  conversion  between  the  susceptible  and  the  persister  is  induced 
by  the  stress  [2],  which  can  be  classified  into  two  categories:  the  stress  from  exogenous  sources,  such 
as  antimicrobial  treatment,  and  the  stress  self-imposed  during  the  bacterial  growth,  such  as  starvation. 
Based  on  this,  we  assume  the  conversion  rate  bsp ,  bps  are  functions  of  concentrations  of  nutrient  and 
antimicrobial  agents: 


bsp  —  bsp(c,  d) ,  bps  —  bps(C)  d). 


(10) 


For  conversion  rate  bsp ,  we  note  that  there  should  be  two  separate  leading  order  terms,  which  rep¬ 
resent  the  effect  of  nutrient  and  antimicrobial  agents,  respectively,  since  these  factors  represent  distinct 
stresses  and  conversion  occurring  at  the  existence  of  the  stress.  For  instance,  nutrient  depletion  alone 
would  induce  persister  formation  [31],  [5],  [6];  without  nutrient  depletion,  persister  formation  would 
also  be  induced  in  response  to  the  antimicrobial  stress  [14].  Thus,  we  propose  bsp  as  follows: 


bsp  — 


®spl 


vspc 


]cnc  rn 
r^spc  \  L' 


+  t>sp2 


(jnd 


Kpd  +  dnd 


1  - 


fap 


bp:  max 


(11) 


where  the  first  term  describes  the  fact  that  bacteria  turn  into  dominant  persisters  due  to  nutrient  depletion 
[31],  the  second  term  depicts  the  ’’stress”  caused  by  antimicrobial  agents  that  promotes  the  conversion 
[14],  and  the  last  term  represents  the  carrying  capacity  for  the  persister.  Here  bsP:  1  and  bsP: 2  represents 
the  maximum  conversion  rates  from  the  susceptible  to  the  persister  due  to  nutrient  depletion  and  the 
antimicrobial  agent  distribution,  respectively.  Although  it  is  plausible  to  propose  one  more  term  in  bsp , 
which  depends  both  on  nutrient  and  antimicrobial  agents,  this  correlation  is  perceived  as  a  high-order 
term  in  this  model;  therefore,  we  omit  it  here  for  simplicity. 

For  the  conversion  rate  bps ,  it  is  assumed  nonzero  only  if  both  the  nutrient  starvation  stress  and 
the  antimicrobial  stress  are  under  certain  threshold  values.  Namely,  it  should  possesses  the  following 
properties:  it  is  zero  when  antimicrobial  agents  level  is  high  enough,  since  a  biofilm  with  persisters  is 
tolerable  to  antimicrobial  treatment,  and  monotonically  increases  to  a  constant  level  as  the  concentration 
of  antimicrobial  agents  drops  below  the  threshold  value,  since  biofilms  are  observed  to  recover  after  an 
antimicrobial  treatment.  In  addition,  the  availability  of  nutrient  can  facilitate  this  process  [6].  One  simple 
rate  is  thus  proposed  as  follows 


bps  —  ^ps,max 


rrnc 


b,md 

psd 


k™sc  +  cm c  k™dd  +  dmd  ’ 


(12) 


where  bpS: max  is  the  maximum  conversion  rate;  kpsc  and  kpsd  are  the  half-saturation  constants,  and  mc 
and  rrid  are  parameters  controlling  the  rate  of  transition  in  the  switch  function. 
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2.5.2  Reactive  kinetics  of  dead  bacteria 


Besides  the  live  bacteria,  the  volume  fraction  of  the  dead  bacteria  is  also  tracked  in  this  model.  It  is 
assumed  the  dead  bacteria  stay  within  the  biofilm  in  the  time  scale  that  this  model  is  valid.  Meanwhile, 
we  assume  some  dead  bacterial  cells  are  attached  to  the  biofilm  acting  as  EPS  while  others  are  converted 
into  solvent  due  to  cell  lysis  [4]  at  certain  specified  rates  in  the  time  scale  of  our  interest, 


9bd  (y bs 


K 2 

^  sd 


Ksd  +  C2 


+ 


C3d 
K3  +  d 


)06s  H-  T'bp&bp  T'dp&bd  ^bd^bdi 


(13) 


where  r ^ ,  rpp  are  the  natural  death  rates  given  before,  C3  is  the  antimicrobial  killing  rate  of  susceptible 
bacteria  with  the  half  saturation  rate  K%,  r^p  represents  the  conversion  rate  from  dead  bacteria  to  EPS, 
and  is  due  to  the  decomposition  of  dead  bacteria  into  solvent. 


2.5.3  Reactive  kinetics  of  EPS 


EPS  is  basically  produced  by  susceptible  bacteria  and  persisters  as  well  as  converted  from  dead 
bacteria.  Over  the  time,  EPS  can  be  dissolved  into  solvent  due  to  reactive  effects.  So,  the  growth  rate  for 
the  EPS  is  proposed  as  follows 


9p  —  [ 


C5C 

k5  +  S 


+ 


C6c 

Kq  +  C 


4>bP](i  -  cp 


^p,  max 


T dp&bd  T'p&pi 


(14) 


where  C5,  Cq  are  the  growth  rate  of  EPS  due  to  susceptible  bacteria  and  persisters,  respectively,  K5,  Kq 
are  half  saturation  constants,  </>p5max  is  the  maximum  volume  fraction  that  EPS  can  achieve  in  the  biofilm, 
and  rp  represents  the  dissolving  rate  of  EPS  into  solvent. 


2.5.4  Reactive  kinetics  of  nutrient 


We  assume  the  nutrient  is  consumed  by  the  live  bacteria  only  and  the  nutrient  decay  rate  is  propor¬ 
tional  to  the  bacterial  volume  fraction: 


9c  —  CV(06s  +  fd^^bp) 


C 

K7  +  c ’ 


(15) 


where  /12  is  a  measure  of  the  nutrient  consumption  efficiency  of  the  persister  cells  relative  to  the  suscep¬ 
tible  cells,  K7  is  a  saturation  constant  and  C7  parametrizes  the  consumption  rate  of  nutrient. 


2.5.5  Reactive  kinetics  of  antimicrobial  agents 

The  antimicrobial  concentration  depends  on  the  live  bacterial  concentration  as  well  as  the  EPS  con¬ 
centration.  It  is  ’’absorbed/consumed”  by  the  live  bacteria  and  possibly  diluted  by  the  EPS  via  chemical 
reactions.  Thus,  the  presence  of  both  bacteria  and  EPS  can  reduce  the  concentration  of  antimicrobial 
agents.  We  thus  propose  the  decay  rate  of  the  antimicrobial  agents  as  follows: 

9d  —  ~{C^(f)bs  +  Cgfibp  +  C\o(f)bd  +  Cn4>p)— — —  +  /d(x,  t)  (16) 

i\8  +  9 

where  is  a  half  saturation  constant,  Cg  is  the  decay  rate  of  antimicrobial  agents  due  to  the  drug- 
susceptible  bacteria  interaction,  Cg  is  the  decay  rate  due  to  drug-persister  cell  interaction,  and  C±o  is  that 
due  to  the  drug  and  dead  bacterial  cell  interaction,  while  C\\  is  that  due  to  drug-EPS  interaction.  Here 
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fd  is  the  source  term  of  antimicrobial  agents,  which  represents  the  injection/supply  rate  of  antimicrobial 
agents  into  the  biofilm,  given  by 


1  (x— xo)2 

/d(x,t)  =  —  M(t)e - A" 

V7 TO 

where  M(t)  is  a  nonnegative,  periodic  function  between  0  and  Mmax  with  periodic  T  and  S  is  a  small 
parameter  limiting  the  range  of  injection  in  space. 

2.6  Boundary  Conditions 

In  this  study,  we  assume  the  biofilm  is  confined  in  a  cubic  domain:  [0  Lx\  x  [0  Ly]  x  [0  Lz] 
where  Lx ,  Ly,  Lz  are  the  length  in  the  x,  y,  z  direction,  respectively.  The  boundary  conditions  for  the 
model  are  proposed  based  on  the  physical  situation  that  we  intend  to  simulate.  In  this  paper,  we  consider 
two  physical  settings  corresponding  to  two  sets  of  boundary  conditions.  One  is  for  a  quiescent  aqueous 
long  channel  and  the  other  is  for  a  short  water  tube  or  flow  cell  with  a  rectangular  cross-section. 

2.6.1  Boundary  conditions  for  the  long  channel 

To  mimic  the  biofilm  development  in  a  long  water  channel,  both  x  and  z  directions  are  assumed 
periodic.  In  the  y  direction,  no-flux  boundary  conditions  are  imposed, 


[cvs(f)s  -  Ds(j)s\7c]  ■  11^=0, Lj, 

=  0, 

[dvsc()s  -  Dd(j)sVd ]  •  n\y=o,Ly 

=  0, 

^ (f>i  •  n|  y=0:Ly 

=  0, 

i  =  6s, 

S  F 

(v</>i  -  A <&V— )  •  n|j/=0;Ly 

0(pi 

=  0, 

i  =  65,  bd,p. 

For  the  average  velocity  we  impose  the  no-slip  condition  on  the  solid  walls  v|y=o,Ly  =  0.  We  can  also 
impose  a  nutrient  feeding  condition  c\y=Ly  =  c*  in  place  of  the  zero-flux  condition  in  certain  parts  of  the 
boundary  as  needed. 


2.6.2  Boundary  conditions  for  the  short  flow-cell 

For  the  case  of  biofilms  in  a  short  water  tube/flow  cell,  periodic  boundary  conditions  are  imposed 
only  in  the  z  direction.  The  x  axis  is  assumed  to  align  with  the  inlet-outlet  direction.  The  inlet  velocity 
at  x  =  0  is  given  by  vq  =  (poyO-  —  y),  0, 0),  where  po  is  a  prescribed  pressure  gradient.  By  assuming 
that  the  solvent  has  already  reached  a  steady  state  while  flowing  out  of  the  cell  at  x  =  Lx,  we  prescribe 
vx  =  0  at  the  outlet  end.  For  the  nutrient  concentration  c,  we  impose  the  feeding  boundary  condition  at 
x  =  0  as  c  =  co(y)  and  cx  =  0  is  assumed  at  x  —  Lx.  For  the  biomass  components,  we  impose  no-flux 
boundary  conditions  in  the  x  direction, 


V c f>i  ■  n\x=o,Lx 

^  f 

)  ’  nU=0 ,LX 

0(Pi 


0,  i  =  bs,bp,bd,p, 
0,  i  =  bs,bp,p,bd. 


We  remark  that  this  boundary  condition  works  for  a  limited  time  frame  before  the  biomass  reaches  the 
boundary.  Beyond  that  moment,  it  must  be  modified.  In  the  y  direction,  we  impose  no-flux  boundary 
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conditions  for  nutrient, solvent,  antimicrobial  agents,  and  biomass  components,  respectively, 


[cvs(f)s  -  Ds(f)s\7c }  ■  n\y=o,Ly 

=  0, 

[dvs(f)s  -  Dd4>sS7d]  ■  n\y=0,Ly 

=  0, 

V(/>*  •  nl^o,^ 

=  0, 

i  =  bs,  bp,  bd,p, 

(v</>i  -  •  n|j/=0;Ly 

0(pi 

=  0, 

i  =  bs,  bp,  bd,p. 

and  the  no-slip  boundary  condition  is  proposed  for  the  velocity  on  the  solid  walls:  v| y=o,Ly  =  0. 


2.7  Nondimensionalization 

Let  to  and  h  represent  the  reference  time  and  length  scale,  respectively.  We  use  these  two  character¬ 
istic  scales  to  nondimensionalize  the  variables  and  equations.  In  this  study,  we  approximate  the  mobility 
parameter  by  a  single  constant  A.  The  nondimensional  variables  are  then  given  by 


Pj 0 


A  =  1?.  r,  =  ^,  r2  = 


_  _  _  /n  _  j  u  r  —  —  rl  =  — 

h  ’  '  poh2^  pQh2  ’  _  co’  do  ’ 


^  1  po/i4  ’  _  2  poh2 

f)  _  Dsto  f)  7  _ 

Us  ~  h?  '  Ud  ~  h2  ’ 
iV 


^kTto  ~  _  ,  p,  ,  p^ 

P  ~  (?SP0  +(PnP0’ 


(17) 


Ci  =  Cito,i  =  1,2,..., 11  -Ki  =  =  2, 5, 6, 7,  spc, psc\  Kj  =  =  3,8 ,spd,psd 


CO 


dp 


hi  —  bito ,  i  =  spl,  sp2,psmax;  T{  =  i  bp,  bd,p. 


where  co,o?o  denote  the  characteristic  substrate  concentration  of  nutrient  and  antimicrobial  agents,  re¬ 
spectively. 

For  simplicity,  we  drop  the  symbol  of  tilde  A  The  nondimensionalized  PDEs  governing  the  biofilm 
system  are  summarized  as  follows: 


p(  f +v-Vv) 

—  V  •  (/Pb^b  “1“  (ppTp  “I-  (ps^i 

V- v 

=  0, 

§-t(t>bs  +  V  •  (06sv) 

V  •  (A (pbs^P  9bs)  9bsi 

Qf<fibp  “I-  ^  '  (^ftpV) 

—  V  •  (A.(pbpW flbp)  “I-  9bpi 

m^bd  +  v  •  (</>Mv) 

=  V  •  (A <pbdS?  9bd)  +  9bd , 

+  v  •  Opv) 

=  V  •  (AcppW /dp)  9p, 

^  +  V-(cvs«/>s) 

=  V  •  (Ds<f>sVc)  +  gc, 

¥  +  V  •  (dvs<ps) 

=  V  •  (De(psV d)  +  gd. 

(18) 
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where  the  reactive  terms  for  each  components  are  given  respectively  as  follows 


9bs  — 

9  bp  — 

9bd 

9p  — 
9c  = 

9d 


C2C 


(i- 


4^bs 


-6, 


+  t>V 


rbsKli  ,  Cad 


f'os^^sd  1 

-  \Jr  2  .  -  + 


i^l+cV"  06a,max^6a+0min  V  i^+c  '  ^3+^ 

^4C  ^ _ A-h  /A,  /A,  _  fV,  -L  C12^ 


)06s> 


K2+c  Va  0fep,max  )^&P  +  bsp4>bs  bps(j)bp  {r^p  +  rf+K3  )06p? 

2  7 

(^6s  x2d+c2  i^3' "h  ^bp^bp  —  T'dp^bd  ~  T'bd&bdi 
[ K5 +c (t)^>s  Ke+c&bplO-  —  Cp(f)p)  +  Tdp4>bd  —  TpCfrp, 
i^bs  +  fd^^bp)  kJ+c  ’ 


—  (C$(/)bs  +  Cgfop  +  Ciofibd  +  Cn^)  K^+d  +  ^=M(t)e  5  ; 


some  of  the  rates  are  prescribed  specifically  as 


bSp  —  [  &sol  7,nc  Tin.  +  6S 


^A^c+^c  ^  sp2k^d+dnd 


j Lmd 
psd 


(1 


(ft  bp 


bps  -  bps, m^krn^+cmc  k^d+dmd- 


(19) 


(20) 


3  Numerical  Methods 


We  will  solve  the  partial  differential  equation  system,  consisting  of  the  governing  equations  of  mo¬ 
mentum,  continuity,  biomass  and  functional  components,  numerically  in  3D  space  and  time  by  a  semi- 
implicit,  finite  difference  method.  We  now  discuss  its  discretization  briefly.  For  simplify,  in  the  follow¬ 
ing,  we  use  the  symbol  of  over-line  for  the  extrapolated  data,  vn+1  =  2vn  —  vn_1. 

Overall,  we  have  10  coupled  PDEs.  We  solve  the  momentum  equation  and  continuity  equation  first 
by  the  Gauge-Uzawa  method  [21].  Recall  that  the  momentum  equation  is  given  by 

dv 

p(—  +  V  •  Vv)  =  V  •  (( hn  +  (f>pTp  +  <j>sTs)  -  [Vp  +  TiV  •  (\/(()n  <g>  V(f)n)]. 

By  adding  a  second  order  term  —  V2v  on  both  side,  it  turns  into 

1  1 

p(-Z  +  V  •  Vv)  -  — V2v  =  -Vp  -  riV2<^nV(/>n  +  V  •  ((j)bTb  +  4>pTp  +  4>sts)  -  —  V2v,  (21) 
Ot  ttea  ^ea 

where  Rea  is  the  averaged  Reynolds  number,  computed  by 


+ 


Pp,avg 

Re„ 


+ 


( fis^avg 


Here  Reb ,  Rep  and  Res  are  the  Reynolds  numbers  for  bacteria,  EPS  and  solvent,  respectively,  and 
06,  p,  s,  avg  are  some  chosen  average  volume  fractions  for  bacteria,  EPS  and  solvent,  respectively.  Then, 
the  Gauge-Uzawa  method  is  adopted  in  three  steps  below.  We  demonstrate  the  numerical  scheme  using 
the  boundary  conditions  for  the  flow  cell  problem. 
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1.  Prediction: 
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2.  Projection: 
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3.  Correction: 
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where 


R' 


72^7n+1V7-7n+1 

+  v  + 


+  V  •  (fb+1rnb+1  +  +p+1++1  + 


fs+1Ts+1) 


=  -i\v2 

Here  s°  =  0  and  v1,  s1, </>£,  </>*  are  computed  by  a  first  order  scheme  and  e  =  0.05  is  used. 

After  the  momentum  and  continuity  equations  are  solved,  the  calculated  velocity  data  are  used  to 
update  the  phase  field  equations.  Since  c/)n  =  q 6^  +  <frbp  +  +  4>p,  we  obtain 


V f^bs  —  ^ l^bp  —  V l^bd  —  V flp  —  riV(V  4>n)  ^2 

If  we  denote 
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the  numerical  schemes  for  the  transport  equations  of  the  effective  biomass  components  are  given  by 
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Analogously,  we  use  the  updated  data  of  the  velocity  and  biomass  to  update  the  concentration  of 
nutrient  and  antimicrobial  agents,  respectively.  The  transport  equations  for  the  nutrient  and  antimicrobial 
agents  are  given  by 


30g+1cn+1-4^cn+0g  1cn  1  vn+ 1  .  V(cn"*"l0n"*"l) 
30s +1^n+1— 40g  rfn+0g  1  dn  1  vn+l  #  ^^n+l^n+l^ 


V  •  (^+1^+1Vcn+1)  +  g™+1, 

V  •  (D™+1^+1Vdn+1)  +  gnd+1, 


where 


9?+1 


e1 


in+ 1\  C7cn+1 


+  ^rbp  )  K7+c™  ’ 

5dn+1  =  -(Csr^+C^+C^+Cn^1)^ 


(29) 


-r^iUrM  JKg+d™- 

The  schemes  are  presented  as  semi-discrete.  The  spatial  discretization  is  carried  out  using  the  central 
difference.  The  boundary  conditions  are  done  using  the  biased  finite  difference  at  the  boundary  coupled 
with  the  central  difference  scheme  for  the  equations.  Grid  refinement  tests  are  carried  out  and  nearly 
second  order  convergence  in  both  time  and  space  are  achieved. 


4  Numerical  Results  and  Discussion 

We  note  that  there  are  many  time  and  length  scales  in  this  hydrodynamic  model.  In  this  paper,  a  few 
selected  regimes  are  studied  which  we  think  are  important  for  the  biofilm  antimicrobial  agent  interaction. 
We  consider  two  characteristic  time  scales  for  to-  the  growth  time  scale  set  at  to  =  103  seconds,  and, 
when  simulating  the  biofilm  treatment  by  antimicrobial  agents  with  hydrodynamic  flows  involved,  the 
antimicrobial  agents  and  biofilm  interaction  time  scale  chosen  as  to  =  10  seconds.  We  remark  that  the 
time  scale  for  biofilm  treatment  can  vary  from  minutes  to  hours,  or  even  days  depending  on  the  choice 
of  the  antimicrobial  agents. 

We  summarize  all  the  model  parameters  used  in  the  current  study  in  table  (1)  with  their  respective 
references  therein.  All  the  parameter  values  used  in  following  simulations  are  chosen  from  table  (1), 
unless  noted  otherwise. 


4.1  Dynamics  of  reactive  kinetics 


Notice  that  the  hydrodynamic  model  is  consisted  of  10  coupled  partial  differential  equations  and  is 
quite  complicated.  Thus,  it’s  advisable  to  study  the  bulk  dynamics  of  the  reactive  kinetics  first,  before 
further  investigating  the  hydrodynamics.  Hence,  by  homogenizing  the  spatial  effects,  our  model  can 
effectively  reduces  to  a  system  of  coupled  ordinary  differential  equations  for  a  spatially  homogeneous 
biofilm  system: 
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We  note  that,  in  our  spatial-temporal  hydrodynamic  model,  nutrient  is  transported  into  the  biofilm  by 
diffusion.  The  last  term  in  nutrient  equation  above  is  modified  to  account  for  this  feature,  where  co  is  the 
nutrient  concentration  in  the  surroundings.  The  reactive  dynamics  of  the  system  of  equations  is  studied 
using  Matlab  ODE  solvers. 

When  the  biofilm  grows  without  antimicrobial  treatment,  all  bacteria  reproduce  provided  that  the 
nutrient  supply  is  sufficient.  Since  the  susceptible  cells  are  more  metabolically  active,  their  population 
grows  exponentially  until  it  reaches  the  carrying  capacity  of  the  environment.  At  the  same  time,  there 
exist  dynamics  of  conversion  between  the  susceptible  and  the  persister  cells.  Here,  we  simulate  the 
spatially  homogeneous  biofilm  development  by  varying  the  maximum  conversion  rates  between  the  sus¬ 
ceptible  and  the  perisister  cells,  i.e.,  bsp i  and  bps^m ax,  respectively.  The  result  is  shown  in  figure  1  and 
figure  2.  Note  that  bsp 2  doesn’t  have  an  effect  in  this  case  since  d  —  0.  As  shown  in  Figure  1,  with  a 
higher  conversion  rate  bsp  1,  the  persister  has  a  higher  volume  fraction  in  steady  state.  A  higher  volume 
fraction  of  the  persister  is  also  observed  by  lowering  conversion  rate  bps  in  Figure  2.  The  conversion  rates 
mainly  affect  the  susceptible  and  persister  populations  and  have  little  impact  on  the  other  components  in 
the  biofilm  as  expected.  This  study  shows  that  the  growth  dynamics  of  the  persister  is  highly  dependent 
on  the  conversion  rates  between  the  susceptible  and  the  persister. 

The  dynamics  of  biofilm  being  treated  by  antimicrobial  agents  are  summarized  in  Figure  3  with 
respect  to  both  constant  dosing  and  periodic  dosing  strategies.  If  we  keep  dosing  constantly,  see  Figure 
3,  the  biofilm  growth  could  be  controlled/contained  but  not  be  eliminated,  as  the  susceptible  cells  could 
be  eradicated  but  the  persisters  are  affected  only  slightly.  However,  if  we  dose  the  sample  periodically, 
that’s  one  dosage  for  a  while  and  ceases  the  dosage  for  a  period  and  then  repeat  the  process  again,  the 
bacteria  in  the  biofilm  can  eventually  be  eradicated  completely  shown  in  Figure  3.  These  predictions 
agree  qualitatively  with  the  literature  [8]  in  that  the  periodic  dosing  strategy  is  viewed  as  an  advisable 
way  to  eliminate  the  biofilm. 

One  important  issue  pertinent  to  the  conversion  of  susceptible  and  persister  cells  is  how  these  two 
bacterial  phenotypes  convert  and  how  fast  they  convert  to  one  another.  While  we  can’t  find  any  timely 
experimental  evidence  to  guide  us  in  this  study  apparently,  we  resort  to  a  systematic  study  using  our 
proposed  model.  Due  to  the  lack  of  refined  studies  in  this  direction,  most  of  the  literature,  such  as  [18], 
regard  bsp  and  bps  as  constants.  Besides,  Cogan  [8]  used  nonconstant  and  nonlinear  conversion  rate 
functions: 


bsp  —  bsp  1 


c  T-  ksp 


bps  —  ^ps,max 


1  +  e 


d—dQ 


(32) 


where  bsp  1  and  bpSj max  are  the  maximum  conversion  rates  and  do  is  the  threshold  for  antimicrobial  agents, 
e  is  a  small  parameter.  However,  we  note  that  these  proposals  are  still  not  quite  satisfactory. 

To  benchmark  the  choices  of  the  conversion  rates  adopted  in  our  model,  we  conduct  several  additional 
numerical  simulations  to  compare  the  constant  conversion  rates  between  the  Cogan’ s  proposal  in  [8]  and 
our  proposal  in  equation  (20).  The  results  are  shown  in  Figure  4. 

In  Figure  4  (a),  for  a  well-grown  biofilm  with  depleted  nutrient  in  the  end,  both  the  susceptible  and 
the  persister  cells  gradually  diminish  in  the  case  of  constant  conversion  rates  or  Cogan’ s  proposal  for 
the  conversion  rates.  These  are  in  direct  conflict  with  the  experimental  observations  where  persisters  are 
in  a  dormant  status  [28]  and  could  survival  in  a  nutrient  depleted  environment  [31]  [5],  as  least  for  a 
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quite  long  period  of  time.  In  comparison,  our  proposed  rates  agree  qualitatively  well  with  the  survival 
persisters  in  the  experiments  in  which  the  persister  decays  slightly. 

Figure  4(b)  depicts  a  scenario  in  which  the  biofilm  is  being  treated  with  antimicrobial  agents.  Note 
that  we  only  plot  the  data  with  the  constant  conversion  rates  and  our  proposal,  as  Cogan’s  conversion 
rates  predict  almost  the  same  as  ours  for  the  susceptible  in  this  case.  For  the  constant  conversion  rates,  all 
bacteria  are  eventually  killed.  However,  this  is  not  the  case  in  the  reality  since  many  experiments  reported 
that  a  small  portion  of  bacteria  would  survive  the  antimicrobial  treatment,  no  matter  how  long  or  how 
strong  the  antimicrobial  doses  are  kept  [28],  which  is  exactly  what  are  observed  using  the  conversion 
rates  in  equation  (20)  proposed  in  this  paper.  So,  our  proposed  conversion  rates  can  handle  both  persister 
survival  due  to  nutrient  depletion  and  antimicrobial  treatment. 

4.2  Heterogeneous  biofilm  development  without  the  antimicrobial  treatment 

Our  main  interest  in  this  paper  is  the  mechanism  of  biofilm  persistence  to  antimicrobial  agents  cou¬ 
pled  with  heterogeneous  biofilm  colonies.  Before  digging  into  the  issue,  we  would  like  to  benchmark 
the  model  for  biofilm  formation  without  antimicrobial  treatment  as  well  as  demonstrate  the  capability  of 
our  proposed  model  for  capturing  details  of  biofilm  formation  in  space  and  time.  In  particular,  how  EPS 
production  can  impact  on  the  growth  of  spatial-temporal  heterogeneous  structures  in  biofilms  and  cell 
starvation  due  to  nutrient  depletion  influences  biofilm  structures  during  biofilm  formation  will  be  studied 
in  detail. 

4.2.1  Heterogeneous  biofilm  structures  regulated  by  EPS  production 

It’s  commonly  perceived  that  EPS  plays  an  essential  role  in  the  formation  and  maintenance  of  the 
structural  stability  of  biofilms  [29]  and  that  EPS  production  rate  could  affect  the  biofilm  structure  dramat¬ 
ically.  In  order  to  study  the  role  of  EPS  production  in  the  process  of  biofilm  formation,  3-D  numerical 
simulations  are  conducted.  Our  numerical  study  agrees  qualitatively  with  the  spatial  variation  in  EPS 
production  [9],  which  leads  to  gradients  in  osmotic  pressure  and  contributes  to  pattern  formation  of 
mushroom  or  tower  shaped. 

Figure  5  depicts  two  numerical  simulations  with  respect  to  varying  EPS  growth  rates.  It  reveals  that 
the  role  of  EPS  production  is  to  expand  the  biofilm  colony  such  that  the  higher  the  EPS  production  rate 
is,  the  thicker  the  biofilm  colony  can  grow.  The  biofilm  growth  without  EPS  production  yields  a  slower 
growth  of  the  biofilm  colony. 

These  numerical  evidence  validates  the  notion  about  the  role  of  EPS  production  in  biofilm  forma¬ 
tion  [37].  Though  higher  EPS  production  would  consume  more  nutrient,  while  in  the  meantime  limits  the 
supply  of  nutrient  to  bacteria.  In  the  case  where  the  nutrient  supply  is  sufficient  by  diffusive  transporta¬ 
tion,  bacteria  would  not  be  affected  negatively.  However,  higher  EPS  production  could  help  to  expand 
the  biofilm  colonies,  as  a  result,  bacteria  could  access  more  nutrient  with  larger  surface  contact  with  the 
solvent,  which  in  turn  help  bacteria  to  grow. 

4.2.2  Effect  of  nutrient  depletion  on  heterogeneous  biofilm  structures 

In  addition  to  EPS  production,  other  factors  may  affect  the  structure  of  biofilm  colonies  during  their 
development  as  well,  one  of  which  is  the  bacterial  death  factor.  There  are  several  causes  that  could 
lead  to  bacterial  death.  These  include  programmed  cell  death  [27],  metabolic  toxins  [30],  competition 
of  different  species  [34]  and  starvation  [7] .  Starvation  is  the  phenomenon  about  bacterial  death  due  to 
nutrient  depletion,  i.e.,  the  bacterial  cells  are  starved  to  death  while  the  nutrient  nearby  is  exhausted. 
At  the  presence  of  cell  death,  it  has  been  observed  that  biofilms  could  form  a  hollow  structure  (void  or 
channel)  due  to  the  decomposition  of  the  dead  bacteria  into  solvent  and  EPS  in  the  biofilm  [24]. 
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We  carry  out  a  numerical  experiment  coupled  with  the  modified  decay  rate  shown  in  Figure  6,  where 
void  structure  is  observed,  as  dead  bacteria  and  parts  of  EPS  are  dissolved  into  the  solvent.  Since  the 
nutrient  is  supplied  on  the  top,  the  level  of  nutrient  is  low  at  the  bottom  part  of  the  domain,  as  shown 
in  figure  7(c).  Initially,  when  nutrient  supply  is  still  sufficient,  the  susceptible  bacteria  grows  homoge¬ 
neously.  Later,  bacteria  at  the  bottom  have  less  access  to  nutrient.  As  it  could  be  observed,  this  results 
in  the  competition  for  nutrient  leading  to  partial  death  of  the  susceptible  bacteria.  Thus,  the  live  suscep¬ 
tible  bacteria  is  mainly  distributed  in  the  upper  part  of  biomass,  where  there  are  more  access  to  nutrient. 
Besides,  the  susceptible  is  more  metabolically  active,  as  could  be  seen  that  the  nutrient  consumption  rate 
is  much  higher  at  the  biomass  surface  in  figure  7(b),  where  it  is  consumed  by  susceptible  bacteria.  In 
the  meanwhile,  the  persister  mainly  distributed  in  the  lower  part,  where  there  is  little  access  to  nutrient, 
since  it  is  in  a  dormant  state  [28],  saying  metabolic  inactive.  The  biomass  flux  is  mainly  at  the  top  sur¬ 
face  of  biofilm  colony,  which  shows  that  nutrient  competition  drives  the  biofilm  to  expand.  Besides,  the 
averaged  velocity  is  much  stronger,  where  there  are  more  volume  of  biomass. 

In  figure  6(i),  we  observe  the  mushroom  shape  of  the  overall-biomass  growth  morphology.  By  adding 
the  feature  of  death  from  starvation  and  competition  for  nutrient,  our  model  could  capture  heterogenous 
biofilm  structure  qualitatively.  We  note  that,  besides  nutrient  competition,  hydrodynamic  stress  would 
also  affect  the  structure  of  biofilm,  as  well  as  the  distributions  of  its  components.  In  this  simulation, 
biofilm  is  grown  in  a  quiescent  environment,  thus  hydrodynamics  stress  is  not  as  effective  as  nutrient 
competition  in  contributing  heterogeneous  biofilm  formation. 

4.3  Biofilm  Development  and  the  impact  of  antimicrobial  agents 

The  concept  that  biofilm  contains  both  persisters  and  susceptible  bacterial  cells  has  been  around  for 
a  while.  Lewis  [28]  has  hypothesized  that  persisters  are  the  main  reason  for  persistent  mechanism  of 
bacteria  cells  to  antimicrobial  agents  in  biofilms.  In  addition  to  the  existence  of  persisters,  some  argues 
that  EPS  acts  as  a  viscous  gel  that  could  prevent  antimicrobial  agents  from  penetrating  deep  into  biofilms 
to  access  the  internal  bacterial  cell;  while  others  propose  that  EPS  could  dilute  the  concentration  of 
antimicrobial  agents  by  reacting  with  it  [42].  Besides,  the  dead  bacteria  may  also  act  as  a  shield  for  the 
biofilm. 

In  the  following  subsections,  we  study  the  mechanism  that  makes  bacteria  living  within  a  biofilm 
more  tolerable  than  planktonic  bacteria  with  respect  to  antimicrobial  treatment.  Later  in  this  section, 
the  effect  of  the  conversion  rate  between  susceptible  and  persister  bacteria  for  heterogeneous  biofilms 
will  be  discussed.  Specifically,  we  conduct  3-D  numerical  simulations  using  our  model  to  reveal  the 
antimicrobial  treatment  of  heterogeneous  biofilms  and  attest  the  hypothesis  of  the  existence  of  persisters 
and  slow  penetration  of  antimicrobial  agents  due  to  EPS,  as  well  as  dead  bacteria.  Then,  we  show  that 
the  periodic  dosing  strategy  is  indeed  a  better  way  compared  with  constant  dosing  in  eradicating  the 
bacteria  in  heterogeneous  biofilms.  For  convenience,  the  initial  biomass  distribution  for  later  simulations 
are  summarized  in  figure  8.  We  note  that  higher  EPS  ratio  in  grown  biofilm  is  usually  observed  in  hy¬ 
drodynamic  environment.  And  we  set  the  volume  fraction  of  persister  bacteria  to  be  zero  in  second  case, 
in  order  to  depict  in  details  the  converting  process  from  the  susceptible  to  persister  during  antimicrobial 
treatment  in  later  simulations. 

4.3.1  Antimicrobial  treatment  and  biofilm  control 

To  develop  a  proper  strategy  for  biofilm  control,  a  clear  understanding  of  the  disinfecting  process 
in  heterogeneous  biofilms  is  necessary.  With  the  3-D  numerical  tool,  we  investigate  the  disinfecting 
process  in  a  long  fluid  channel  as  well  as  a  short  flow  cell,  during  which  various  dosing  positions  will  be 
examined. 
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As  shown  in  Figure  9,  we  mimic  the  situation  where  the  biofilm  develops  in  a  tank  with  quiescent 
solvent  and  antimicrobial  agents  supplied  in  a  well-mixed  form  with  the  pure  solvent  from  the  top  bound¬ 
ary.  In  reality,  this  scenario  agrees  with  the  case  of  spraying  antimicrobial  solution  on  the  biofilm  surface. 
Given  a  well-grown  biofilm,  once  we  start  dosing  antimicrobial  agents,  the  susceptible  cells  will  be  killed 
drastically.  However,  the  volume  fraction  of  the  persisters  retains  at  a  stable  level,  as  well  as  EPS.  In 
other  words,  this  dosing  method  fails  to  eradicate  the  grown  biofilm  colony  completely,  which  is  always 
the  case  in  reality. 

The  second  numerical  experiment  mimics  the  situation  of  disinfection  in  a  flow-cell,  where  the  an¬ 
timicrobial  agents  are  well-mixed  with  the  solvent  and  flown  in  from  the  inlet  boundary.  In  reality,  this 
approximates  the  disinfection  of  biofilms,  growing  in  a  river  bed  or  in  a  water  pipe.  The  simulation 
results  are  shown  in  Figure  10.  Unlike  the  case  alluded  to  earlier,  the  distribution  of  dead  bacteria  is 
heterogeneous,  namely,  more  susceptible  bacteria  are  killed  facing  the  inlet  boundary  than  facing  the 
downstream  due  to  the  heterogeneous  distribution  of  antimicrobial  agents.  In  the  meanwhile,  some  sus¬ 
ceptible  bacteria  are  converted  into  persister  bacteria,  as  the  existence  of  antimicrobial  agents  triggered 
this  conversion.  In  order  to  depict  more  details,  2D  slices  at  z  =  0.5  are  shown  in  figure  1 1.  As  could  be 
observed,  antimicrobial  agents  are  distributed  heterogeneously,  as  well  as  its  consumption  rates.  Due  to 
the  limitation  of  antimicrobial  diffusive  rate  into  biofilm,  they  are  not  consumed  totally  at  the  inlet  bound¬ 
ary,  but  convected  through  the  domain.  Effective  hydrodynamic  pressure  is  high  at  the  inlet  boundary  and 
low  in  the  channels  of  biofilm  colonies.  We  note  that  there  is  biomass  flux  towards  the  outlet  boundary, 
as  shown  in  figure  11(f).  However,  the  biofilm  is  in  a  laminar  flow,  thus  its  the  disinfection  process  that 
dominate  the  dynamics.  Similar  with  the  previous  case,  persisters  can  survive  this  treatment,  in  other 
words,  this  dosing  method  also  fails  to  eradicate  the  biofilm  colony  satisfactorily. 

We  note  that,  to  the  best  of  our  knowledge,  most  mathematical  models  available  treating  antimicro¬ 
bial  agents  as  well-mixed  with  solvent  and  thereby  describe  the  antimicrobial  agents  using  a  spatially 
homogeneous  differential  equation  (ODE).  In  real  applications  however,  we  must  take  into  account  the 
heterogeneous  distribution  of  antimicrobial  agents  and  biofilm  structures.  Therefore,  we  consider  an  in¬ 
jection  of  antimicrobial  agents  into  the  solvent  by  a  needle.  After  the  injection,  the  antimicrobial  distribu¬ 
tion  within  the  biofilm  ambient  fluid  ensemble  is  dominated  by  convection  and  diffusion  of  antimicrobial 
agents  through  various  media.  As  far  as  we  know,  it  is  the  first  mathematical  biofilm  model  and  simula¬ 
tion  in  the  literature  to  show  this  full  3D  view  of  disinfection  process  by  antimicrobial  injection.  Figure 
12  depicts  the  simulation  details.  During  the  process,  persisters  are  converted  from  susceptible  bacteria 
and  the  killing  of  bacteria  is  asymmetric  and  heterogeneous  in  space  and  time.  The  concentration  of  the 
dead  bacteria  close  to  the  source  of  the  antimicrobial  agents  is  apparently  much  higher  than  that  on  the 
other  side  of  the  biofilm  colony.  The  persister  cell  concentration  doesn’t  drop  at  the  dosing  position,  in 
other  words,  persisters  survive  this  treatment  as  well. 

Despite  that  constant  dosing  shown  in  the  simulations  can  not  eradicate  the  bacteria  in  the  biofilm 
colony  completely,  our  model  suggests  that  targeted  disinfection,  i.e.,  disinfection  by  injection  to  the 
specific  position,  is  more  efficient  and  environmental-friendly  since  it  requires  less  amount  of  antimicro¬ 
bial  agents  to  achieve  the  same  effect,  especially  for  biofilms  grown  in  an  aqueous  system,  where  the 
flow  velocity  is  not  too  small  to  be  ignored.  With  the  development  of  nanotechnology,  targeted  delivery 
of  antimicrobial  agents  directly  to  the  biofilm  is  becoming  an  reality  and  is  shown  to  be  more  efficient. 
Moreover,  since  EPS  could  prevent  antimicrobial  agents  from  diffusing  away  from  the  biofilm,  the  local 
antimicrobial  concentration  could  be  maintained  for  longer  period  of  time. 

We  note  that  the  diffusion  rate  of  antimicrobial  agents  is  partially  responsible  for  the  heterogeneous 
distribution  of  dead  bacteria.  If  the  diffusion  rate  of  antimicrobial  agents  is  high,  the  distribution  of  dead 
bacteria  would  be  homogeneous  with  respect  to  biofilm  colonies.  But  If  the  diffusion  rate  is  small,  in 
other  words,  the  transport  in  the  biofilm- solvent  system  is  convection  dominated,  the  density  of  dead 
bacteria  would  be  highly  heterogeneous  with  higher  concentration  near  the  dosing  position. 
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Observed  from  the  disinfection  process  above,  every  constant  dosing  case  fails  to  eradicate  all  the 
bacteria  in  one  long  dose.  Therefore,  a  straightforward  question  is:  what  happened  to  the  biofilm  after 
the  cessation  of  an  antimicrobial  dosing?  One  simulation  result  is  shown  in  Figure  13  with  the  initial 
biofilm  colony  from  the  experiment  in  Figure  9.  After  dosing  stops  and  the  concentration  of  antimicrobial 
agents  drops  below  the  threshold,  persisters  become  active  from  dormant  status  and  convert  back  into 
susceptible  cells  for  proliferation.  This  leads  to  the  relapse  of  biofilm  growth  after  certain  period  of 
time.  The  newly  grown  biofilm  is  morphologically  similar  with  previous  colony,  as  could  be  observed. 
Besides,  we  also  absolve  that  the  susceptible  bacteria  mainly  distributed  at  the  surface  of  biomass,  where 
they  have  more  access  to  nutrient.  But,  the  persister  bacteria  is  less  metabolical  active  [28],  staying  at 
the  bottom  of  the  biomass,  where  there  is  less  nutrient. 

4.3.2  Slow  penetration  effects  due  to  EPS  and  dead  bacteria 

It’s  commonly  known  that  bacteria  living  in  a  biofilm  are  more  persistent  to  antimicrobial  treatment 
than  planktonic  bacteria.  Besides  the  existence  of  persister  cells,  reason  for  this  phenomenon  is  that  EPS 
can  acts  as  a  shield  for  trapped  bacteria  and  thereby  could  retard  or  even  prevent  the  antimicrobial  agents 
from  penetrating  deeper  into  the  biofilm  colony.  Another  potential  shield  is  that  some  dead  bacteria 
would  gather  on  some  internal  surface  of  the  biofilm,  which  also  makes  it  harder  for  the  antimicrobial 
agents  to  diffuse  into  the  biofilm  . 

We  investigate  the  phenomenon  of  antimicrobial  penetration  using  our  model  in  3-D  simulations  for 
varying  diffusion  rate  of  antimicrobial  agents  proposed  in  equations  (5), (6), (7).  Initially  we  begin  with 
a  well-grown  biofilm  in  the  center  of  the  domain,  the  simulation  results  for  the  three  cases  with  distinct 
diffusion  coefficients  are  shown  in  Figure  14.  Slow  penetration  effects  due  to  EPS,  as  well  as  dead 
bacteria,  are  observed  during  the  disinfection  process. 

The  total  volume  fraction  of  each  component  after  antimicrobial  treatment  is  shown  in  Figure  15. 
The  difference  is  the  varying  decay  rate  of  susceptible  bacteria  due  to  slow  penetrating  effect  induced 
by  EPS  and  dead  bacteria.  Details  of  the  susceptible  decaying  rate  is  shown  in  Figure  15(A).  Due  to 
the  slower  penetration  of  antimicrobial  agents  at  the  presence  of  EPS  and  bacteria,  the  survival  curve 
drops  slower  than  the  controlled  set,  i.e.  constant  diffusion  rate  proposed  in  equation  (5),  since  more 
susceptible  bacteria  survive  and  they  could  mutate  into  persisters.  This  gives  a  convincing  explanation 
that  EPS  acts  as  a  shield  for  the  bacteria  [16]. 

However,  apparently  from  (E)  in  Figure  15,  each  case  fails  to  eradicate  the  bacteria  in  the  biofilm, 
as  a  small  portion  of  persisters  is  hardly  affected  by  antimicrobial  agents  shown  in  Figure  15(C).  On 
the  other  hand,  our  numerical  results  demonstrate  that  the  slow  penetration  caused  by  EPS  and  bacteria 
distribution  does  slow  down  the  antimicrobial  treatment  of  the  biofilm.  However,  it  only  shows  as  a  factor 
for  retarding  the  treatment  of  biofilm,  but  not  the  essential  component  for  antimicrobial  persistence  of 
biofilm. 

4.3.3  Dosing  strategy 

As  reported  in  [28]  experimentally  as  well  as  predicted  by  our  biofilm  model,  antimicrobial  agents 
fails  to  eradicate  all  the  bacteria  in  one  long  dose.  After  the  dosing  is  stopped,  the  bacterial  growth  can 
relapse.  Thus,  searching  for  an  optimal  dosing  strategy  is  necessary. 

When  dosing  strategies  are  considered,  basically,  there  are  two  ways  to  carry  it  out.  One  is  to  give 
high  concentrated  doses  of  antimicrobial  agents  in  short  time  period  and  the  other  is  to  give  low  concen¬ 
tration  doses  of  antimicrobial  agents  continuously  at  a  long  time  period.  The  question  here  to  discuss 
is  which  way  is  better.  Although  several  research  articles  have  discussed  this  issue  for  homogeneous 
biofilms  [39],  [19],  [20]  and  [8],  it  is  worthy  of  exploring  the  solution  in  a  heterogeneous  biofilm.  In 
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addition,  the  poorly  understood  mechanism  of  biofilm  resistance  to  antimicrobial  agents,  more  practi¬ 
cal  issues  arise  recently,  such  as  the  cost  of  antimicrobial  agents  and  the  requirement  of  environmen¬ 
tally  friendly  antimicrobial  agents  [19],  which  demands  less  usage  of  antimicrobial  agents  to  control  the 
biofilm. 

As  previously  studied,  constant  dosing  could  control  the  bacteria  to  within  a  low  level,  but  it  can 
leave  a  small  portion  of  them  alive,  most  of  which  are  persisters.  Since  the  conversion  rate  of  persisters  to 
susceptible  cells  with  existing  antimicrobial  agents  is  relatively  small,  instead  of  being  killed  by  changing 
into  susceptible  bacteria,  persisters  could  keep  its  number  stable  for  a  quite  long  period  of  time.  Thus, 
after  the  susceptible  cells  are  eradicated,  a  short  period  of  ceased  dosing  is  in  favor  for  facilitating  the 
persister  conversion  into  susceptible  cells,  which  can  then  be  killed  easily  by  the  subsequent  dosages. 
In  Figure  16,  for  a  given,  well-grown  biofilm,  we  dose  antimicrobial  agents  for  10  units  of  time,  where 
to  —  1000  seconds  and  then  cease  dosing  for  20  units  of  time,  and  then  repeat  this  process.  Following 
this  strategy,  both  the  susceptible  and  the  persister  cells  can  be  killed  eventually. 

With  this  hydrodynamic  model,  we  confirm  the  hypothesis  that  periodic  dosing  strategies  give  higher 
disinfection  rate  and  could  eradicate  all  the  bacteria  or  to  reduce  its  concentration  to  under  a  prescribed, 
lower  level  within  a  short  period  of  time  provided  that  the  dosing  period  and  dosing  strength  is  well 
controlled. 

5  Conclusion 

In  this  paper,  we  extend  our  previous  phase  field  models  for  biofilms  [47],  [48]  to  multiphasic  biofilms 
by  taking  into  account  multiple  types  of  bacteria  and  the  effect  of  antimicrobial  agents.  Then,  a  series 
of  three  dimensional  numerical  simulations  are  carried  out  to  investigate  biofilm  development  both  with 
and  without  antimicrobial  treatment. 

For  naturally  growing  biofilms  without  antimicrobial  treatment,  our  model  verifies  the  hypothesis  that 
the  production  of  extra  polymeric  substances  promotes  the  spreading  of  biofilms  by  generating  osmotic 
pressure  [37].  In  addition,  the  void  structure  of  the  biofilm  morphology  due  to  nutrient  depletion  can  also 
be  captured,  which  qualitatively  agrees  with  the  experiment  observation  [15]. 

For  the  antimicrobial  treated  biofilms,  both  antimicrobial  treatment  in  an  infinite-long  channel  and  a 
finite-long  tube,  along  with  varying  dosing  positions  and  strategies  are  investigated.  As  suggested  from 
the  numerical  simulations,  the  mechanism  of  conversion  between  the  susceptible  bacteria  and  persisters 
is  essential  for  the  growth  dynamics  of  persisters. 

Besides  persisters,  EPS  and  dead  bacteria  can  form  obstacles  to  hinder  the  diffusion  of  antimicrobial 
agents,  which,  in  turn,  protects  biofilm  from  being  disinfected  quickly.  As  a  result,  it  provides  more 
time  for  susceptible  bacteria  to  transform  into  persisters.  The  mechanisms  of  retarding  the  antimicrobial 
diffusion  by  EPS  and  dead  bacteria  make  biofilms  even  harder  to  be  disinfected. 

It  is  difficult  to  eradicate  biofilms  completely  by  conventional  means.  Thus,  proper  strategies  are  in 
need  for  biofilm  control.  As  is  observed  from  the  numerical  simulations,  dosing  by  injection  is  much 
more  efficient,  as  well  as  environment-friendly,  than  using  a  nebulizer,  which  delivers  the  antimicro¬ 
bial  agent  to  the  surface  of  the  biomass-solvent  mixture.  Especially,  antimicrobial  agent-carrying  nano 
spheres,  which  could  be  imbedded  within  the  biofilm,  would  be  much  more  efficient,  as  EPS  would  pre¬ 
vent  the  antimicrobial  agents  from  diffusing  away  from  the  biofilm.  Besides  the  dosing  position,  dosing 
strategies  are  also  vital  for  effective  disinfection.  Our  numerical  simulations,  consistent  with  experiments 
reported  in  literature,  confirm  that  the  periodic  dosing  would  make  a  better  biofilm  control  than  dosing 
for  one  time  or  continuous  dosing  [19]. 

In  general,  the  three  dimensional  heterogeneous  biofilm  model  is  capable  of  capturing  details  of 
biofilm  development  and  morphological  changes  during  antimicrobial  treatment.  It  provides  an  in-silico 
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tool  for  analysing  the  mechanism  of  biofilm  persistence  to  antimicrobial  agents  and  deriving  potentially 
optimal  dosing  strategies  for  biofilm  control  or  disease  treatment. 
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Symbol 

Description 

value 

Unit 

Source 

T 

absolute  temperature 

303 

Kelvin 

[47] 

ks 

Boltzmann  constant 

1.38065x  10-23 

m?kgs~ 2 

[47] 

Pn 

biomass  density 

1  x  103 

kg  m~3 

[47] 

Ps 

water  density 

1  x  103 

kg  m~3 

[47] 

h 

characteristic  length  scale 

1  X  10“3 

m 

[32] 

to 

characteristic  time  scale 

10  or  1000 

s 

[47] 

Lx,  Ly  ,  Lz 

size  of  computational  domain 

1  -  2  x  10-3 

m 

[47] 

Pb 

dynamic  viscosity  of  bacteria 

2.7  x  102 

kg  m~1s~ 1 

[25] 

Pp 

dynamic  viscosity  of  EPS 

2.7  x  102 

kg  m~1s~ 1 

[25] 

Ps 

dynamic  viscosity  of  solvent 

1.002  x  10"3 

kg  m~1s~ 1 

[25] 

co 

characteristic  oxygen  concentration 

8.24  x  10“3 

kg  m~3 

[32] 

do 

characteristic  antimicrobial  concentration 

8.24  x  10“3 

kg  m~3 

71 

distortional  energy  coefficient 

8  x  106 

m~ 1 

72 

mixing  free  energy  coefficient 

3  x  1017 

m~3 

X 

Flory-Huggins  parameter 

0.55 

[47] 

A 

mobility  parameter 

1  x  10-9 

kg~1m3s 

N 

generalized  polymerization  parameter 

1  x  103 

[47] 

Dc 

oxygen  diffusion  coefficient 

2.3  x  10“9 

m~2s~ 1 

[41]  ,  [32] 

Dd 

antimicrobial  diffusion  coefficient 

2.3  x  10“8 

m~2s~1 

[41] 

Dpr 

Hinson  model  parameter 

0.02 

[22] 

06s, max 

carrying  capacity  for  susceptible  bacteria 

0.2 

06p,  max 

carrying  capacity  for  persister 

0.02 

0p,max 

carrying  capacity  for  EPS 

0.2 

c2 

susceptible  bacteria  growth  rate 

4  X  10“4 

s-1 

c3 

susceptible  bacteria  decaying  rate 

4  x  10-2 

s-1 

c4 

persist  bacteria  growth  rate 

4  x  lO"7 

s-1 

Ci2 

persister  bacteria  decay  rate 

4  x  10“5 

s_1 

rbs,  max 

flush-out  rate  for  susceptible 

4  x  10“7 

S_1 

nP 

flush-out  rate  for  persister 

1  x  10“7 

s_1 

rbd 

flush-out  rate  for  dead  bacteria 

1.0  x  10“7 

s_1 

t>spl 

transfer  rate 

1  x  lO"5 

s_1 

t>sp2 

transfer  rate 

1  x  10“3 

s  1 

kSpd ,  kSpc 

Monond  constant  for  bsp 

3.5  x  10“4 

kg  m~3 

bps,  max 

transfer  rate  from  to 

4  x  10“5 

s_1 

kpsc ? kpSd 

Monod  constant  for  bps 

3.5  x  10“4 

kg  m~3 

mc,  rrid ,  nc ,  nd,  n 

Hill  parameter 

2 

EPS  growth  rate  due  to  susceptible 

4  x  10“4 

s  1 

c6 

EPS  growth  rate  due  to  persisters 

4  x  lO"6 

s_1 

Tp 

flush-out  rate  for  EPS 

4  x  10“7 

s~1 

rdp 

converting  rate  from  death  bacteria  to  EPS 

4  x  10“7 

s-1 

c7 

nutrient  consumption  rate 

4  x  lO"2 

kg  m~3  s~ 1 

P2 

nutrient  consumption  raters  v.s.  4>bp) 

0.1 

kcd 

Monond  constant  for  r^d 

3.5  x  10“4 

CS,C9 

antimicrobial  agents  consumption  rate 

4  x  10-2 

kg  m~3  s_1 

C\0,  Cn 

antimicrobial  agents  consumption  rate 

4  x  lO"4 

kg  m~3  s~ 1 

K2,K3,K5,Kg,K8 

monod  constant  for  EPS  growth 

3.5  x  10“4 

kg  m~3 

[32] 

Table  1:  The  table  for  parameter  values.  All  parameters  with  resources  are  marked.  The 
others  are  fitted  from  our  previous  work  and  experiment  results. 
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Figure  1:  Growth  dynamics  of  spatially  homogeneous  biofilms.  This  figure  shows  biofilm 
growth  dynamics  with  respect  to  four  selected  conversion  rate  bsp i:  0,  2xl0-5, 2xl0-4, 2x 
10-3.  The  volume  fractions  of  the  susceptible,  persister,  dead  bacteria,  EPS  and  the  concen- 
tration  of  nutrient  and  total  volume  fraction  of  live  bacteria  with  time  are  shown  in  (a)-(f) 
respectively. 
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Figure  2:  Growth  dynamics  of  spatially  homogeneous  biofilms.  This  figure  shows  biofilm 
growth  dynamics  with  respect  to  four  selected  conversion  rates  6pS}max-  4  x  10-5,4  x 
10-4, 4  x  10-3, 4  x  10-2.  The  volume  fractions  of  the  susceptible,  persister,  dead  bacteria, 
EPS  and  the  concentration  of  nutrient  and  total  volume  fraction  of  live  bacteria  with  time 
are  shown  in  (a)-(d)  respectively. 
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Figure  3:  Dynamics  of  spatially  homogeneous  biofilms  being  treated  by  antimicrobial 
agents.  This  figure  shows  two  simulations:  one  with  antimicrobial  agents  dosed  constantly 
and  the  other  with  antimicrobial  agents  dosed  periodically.  The  former  fails  to  kill  the  bac¬ 
teria  in  the  biofilm  while  the  latter  can  eradicate  the  bacteria  in  the  biofilm.  Here  we  suppose 
nutrient  is  sufficient.  The  volume  fractions  of  the  susceptible,  persister,  dead  bacteria,  EPS, 
the  concentration  of  antimicrobial  agents  and  the  volume  fraction  of  live  bacteria  with  time 
are  shown  in  (a)-(f),  respectively. 
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(a)  Biofilm  growth  dynamics  with  different  conversion  rates  bsp,  bps. 


(b)  Biofilm  dynamics  while  treated  by  antimicrobial  agents  with  different  conversion  rates  bsp,  bps 


Figure  4:  A  Comparison  of  biofilm  growth  and  dynamics  after  being  treated  by  antimi¬ 
crobial  agents  with  different  conversion  rates  bsp  and  bps.  In  this  figure,  three  conversion 
model  for  bsp  and  bps  are  used:  the  constant  conversion  rates  bsp  =  bsp i,  bps  =  &ps?max> 
the  conversion  rates  of  equation  (32),  where  do  =  0.1,  e  =  0.01  and  the  conversion  rates 
of  equation  (20),  which  are  represented  in  the  legend  by  ‘Constant’,  ‘Cogan’  and  ‘Mine’, 
respectively,  (a)  A  comparison  biofilm  growth  dynamics  when  nutrient  supply  is  deficient, 
(b)  A  comparison  of  dynamics  for  biofilm  being  treated  by  antimicrobial  agents. 
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(a)  4>hs  at  t  =  0 
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0.012 


0.0090 


0.0060 
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(f)  4>p  at  t  =  212 


J  0.014 
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0.0082 
0.0055 
0.0027 


h  0.040 
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(g)  (j)hs  at  t  =  212 


(h)  4>bP  at  £  =  212 


(i)  (\)pdXt  —  212 


Figure  5:  Biofilm  development  with  varying  EPS  production  rates.  It  shows  that  enhanced 
EPS  production  can  expand  the  biofilm,  (a)-(c)  initial  volume  fractions  of  the  susceptible, 
the  persister  and  EPS;  At  time  t  —  212,  volume  fractions  of  the  susceptible,  the  persister  and 
EPS,  when  the  EPS  production  rate  C5  =  8  x  10-4,  are  shown  in  (d)-(f)  respectively;  (g)- 
(i)  are  the  volume  fractions  of  the  susceptible,  the  persister  and  EPS  without  EPS  prodution 

(C5  =  °)- 
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(g)  (j)p  at  t  =  100 


(h)  fad  at  t  =  100 


(i)  at  t  =  100 


Figure  6:  Void  structures  in  biofilm  colonies  due  to  the  death  of  bacteria  as  a  result  of 
depleted  nutrient.  This  figure  shows  heterogenous  biofilm  growth  dynamics  of  a  single 
biofilm  bump.  In  this  case,  nutrient  is  supplied  on  the  top  and  initially  there  is  a  single 
bump  containing  only  susceptible  bacteria,  where  we  choose  </>p5max  =  0.1,  r m  =  10-3, 
rbs  =  rp  =  10-5.  In  order  to  depict  it  clearly,  all  plots  are  2D  slices  at  z  =  0.5.  (a) 
initial  biofilm  colony;  (b)-(d)  volume  fraction  of  susceptible  bacteria  at  t  =  30,  60, 100 
respectively;  (e)-(f)  volume  fraction  of  persister  bacteria  at  t  =  60, 100  respectively;  (g)-(i) 
volume  fraction  of  EPS,  dead  bacteria  and  total  biomass  at  t  =  100  correspondingly. 
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Figure  7:  Void  structures  in  biofilm  colonies  due  to  the  death  of  bacteria  as  a  result  of  de¬ 
pleted  nutrient.  This  figure  shows  other  2D  slice  plots  (z  =  0.5)  for  simulations  in  figure 
6  at  time  t  =  200.  (a)  hydrodynamic  pressure;  (b)  nutrient  consumeption  rate;  (c)  nutrient 
distribution;  (d)  maximum  shear  stress  due  to  biomass;  (e)  effective  biomass  flux;  (f)  vol¬ 
ume  averaged  velocity;  (g)external  surface  force  (force  derived  by  least  action  principle  in 
moment  equation). 
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(a)  4>n  at  t  =  0 


(b)  4>n  at  t  =  0 


Figure  8:  Initial  Condition  for  later  biofilm  treatment  simulations.  This  figures  show  well- 
grown  biofilms,  which  will  be  used  as  initial  conditions  for  later  simulations:  (a)  volume 
fraction  of  total  biomass  for  a  grown  biofilm  in  an  infinite-long  channel,  where  we  simply 
suppose  the  susceptible,  persister  and  EPS  and  dead  bacteria  has  same  distribution  with 
volume  fraction  ratio  as  3  :  1  :  1  :  0;  (b)  grown  biofilm  in  a  water  tube,  where  the  corre¬ 
sponding  ratios  are  proposed  as  1  :  1.2  :  0  :  0. 
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Figure  9:  Antimicrobial  treatment  process  in  an  infinite-long  channel  using  well-mixed 
antimicrobial  agents  dosed  from  the  top  boundary.  This  figure  shows  susceptible  bacteria 
are  dramatically  eradicated,  leaving  only  persister  cells:  (a)-(c)  susceptible  bacteria  volume 
fraction  at  time  t  =  10,  36, 134  respectively;  (g)-(e)  three  slices  of  susceptible  bacteria  at 
t  =  0, 10  correspondingly;  (f)  antimicrobial  agent  concentration  at  t  =  134  (g)-(i)  dead 
bacteria  and  persister  and  EPS  at  t  —  134,  respectively. 
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(a)  4>hs  at  t  =  2 


(b)  (j)hs  at  t  =  14 


(c)  0bs  at  t  =  24 


(d)  4>bp  sit  =  2 


(e)  fop  at  t  =  24 


(f)  4>bd  at  t  =  24 


Figure  10:  Biofilm  treatment  in  a  finite-long  tube  using  well-mixed  antimicrobial  agents 
dosed  from  inlet  boundary.  Here  the  reference  time  is  to  —  10  seconds  and  the  initial 
biofilm  colony  is  shown  in  figure  8(b).  (a)-(c)  the  volume  fraction  of  susceptible  bacteria  at 
time  t  =  2, 14, 24,  respectively;  (d)-(e)  the  volume  fraction  of  persister  bacteria  at  t  =  2, 24 
respectively;  (f)  the  volume  fraction  of  dead  bacteria  at  t  =  24. 
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X-Axis 
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(c)  Concentration  of  Antimicrobial  Agent  at  t  =  2  (d)  Effective  Hydrodynamic  Presure  at  t  =  2 


(e)  v  at  t  =  2  (f)  0nvn  at  t  =  2 

Figure  11:  Biofilm  treatment  in  a  finite-long  tube  using  well-mixed  antimicrobial  agents 
dosed  from  inlet  boundary.  This  figure  shows  2D  slice  (z  =  0.5)  for  the  simulation  in  figure 
10  at  time  t  =  2.  (a)antimicrobial  agent  consumption  rate;  (b)  total  biomass  distribution; 
(c)  antimicrobial  agent  distribution;  (d)  hydrodynamic  press;  (e)  volume-averaged  velocity; 
(f)  biomass  flux. 
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(c)  (f)bd  at  t  =  1463 
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(d)  datt  =  310 


(e)  4>bs  at  t  =  1463 


(f)  fop  at  t  —  1463 


(g)  at  t  —  1463 


Figure  12:  Biofilm  treatment  by  antimicrobial  agents  injected  through  an  injection  needle. 
A  highly  heterogeneous  distribution  of  dead  bacteria  is  observed  due  to  the  pinpointed  dos¬ 
ing  strategy.  The  initial  profiles  is  given  in  figure  8(a)  and  the  injection  position  is  chosen 
as  (0.6,  0.1,  0.4).  (a)-(c)  provides  the  profiles  of  dead  bacteria  at  varying  time:  t  =  124, 
310,  1463,  respectively;  (d)  2D  slice  (z  =  0.4)  of  concentration  of  antimicrobial  agents  at 
t  =  310;  the  volume  fractions  of  the  susceptible,  the  persister  and  EPS  at  time  t  =  1463  are 
given  in  (e)-(g),  respectively. 
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Figure  13:  Biofilm  relapse  after  the  cessation  of  an  antimicrobial  agents  dosing.  The  biofilm 
is  observed  to  relapse  as  susceptible  bacteria  regrow  into  colonies  due  to  the  persister- 
susceptible  conversion  after  the  treatment.  The  initial  biofilm  colony  is  an  antimicrobial- 
treated  one  taken  from  the  numerical  experiment  in  Figure  9,  where  all  the  susceptible  bac¬ 
teria  are  eradicated,  (a)-(c)  shows  susceptible  bacteria  at  time  t  —  76, 190, 235  respectively; 
(d)-(f)  shows  the  volume  fractions  of  dead  bacteria,  the  persister  and  EPS  at  t  =  235  corre¬ 
spondingly;  (g)-(i)  shows  three  slices  of  the  volume  fractions  for  the  susceptible,  persister 
bacteria  and  EPS  at  t  =  235  respectively. 
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(a)  (j)bs  at  t  =  10 


(d)  (j)bs  at  £  =  30 


(e)  fas  at  t  —  30 


(f)  4>bs  at  t  =  30 


Figure  14:  A  comparison  of  varying  antimicrobial  diffusion  rates  in  a  disinfection  process. 
This  figure  shows  antibiotic  treatment  of  biofilms  in  an  infinite-long  channel  with  antimicro¬ 
bial  agents  dosed  on  the  upper  boundary,  using  three  different  diffusion  rates:  De i,  De 2  and 
De 3.  The  initial  condition  is  shown  in  Figure  8(a).  The  volume  fraction  of  the  susceptible 
bacteria  for  penetrating  rates  Dei,De2  and  De 3  at  t  =  10  are  shown  in  (a)-(c)  respectively; 
and  (d)-(f)  shows  the  corresponding  volume  fractions  of  the  susceptible  bacteria  at  t  —  30. 
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Figure  15:  A  comparison  of  varying  antimicrobial  diffusion  rates  in  a  disinfection  process. 
This  figure  shows  total  volume  fractions  of  each  components  for  simulations  in  figure  14: 
(a)-(d)  total  volume  fractions  of  susceptible,  dead,  persister  bacteria  and  EPS,  correspond¬ 
ingly;  (e)  total  volume  fraction  of  live  bacteria;  (f)  total  concentration  of  antimicrobial 
agents. 
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Figure  16:  Periodic  dosing  strategies  succeed  in  treating  biofilms.  Given  a  well-grown 
biofilm  from  figure  8(a),  antimicrobial  agents  are  dosed  periodically.  Then,  the  live  bacteria 
are  gradually  eradicated.  The  volume  fraction  of  persister  bacteria  at  varying  times:  t  = 
212, 1098, 1782  are  shown  in  (a)-(c)  respectively;  (d)-(f)  shows  the  volume  fractions  for  the 
susceptible,  dead  bacteria  and  EPS,  correspondingly. 
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